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THESIS DISCLAIMER 



The reader is cautioned that computer programs developed 
in this research may not have been exercised for all cases of 
interest. While every effort has been made, within the time 
available, to ensure that the programs are free of 
computational and logic errors, they cannot be considered 
validated. Any application of these programs without 
additional verification is at the risk of the user. 
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I. 



INTRODUCTION 



A. BACKGROUND 

In 1966 W.R. Blischke and A.H. Halpin [Ref. 1] first 
described methods for the computation of confidence intervals 
for circular error probability (CEP) based on first order 
variance estimates. When these procedures were implemented at 
the Johns Hopkins University Applied Physics Laboratory (JHU- 
APL) it was found that under certain conditions the resulting 
confidence intervals for CEP were smaller than expected. As 
a result R.C. Ferguson, K.V. Kitzman and P.B. Jackson [Ref. 2] 
have extended the Blischke-Halpin procedures to create a 
second order variance estimate. Testing conducted by Kitzman 
[Ref. 3] has confirmed that this method produces a more 
accurate estimate for the variance of CEP. 

Our goal will be to extend this procedure to the 3- 
dimensional case to obtain a second order estimate for 
variance of spherical error probability (SEP) . 

B. METHODOLOGY 

We begin by assuming that missile detonations are 
distributed as trivariate Gaussian. In the past, analysis of 
ballistic missile test firing data has failed to disprove this 
assumption. If the mean and variance are given by: 



1 



n = 






s = 



°11 °12 °13 

a 21 °22 °23 

a 31 °32 °33 



Then the probability of a detonation occurring within a 
distance R of the origin is given by: 

p, 2) = ff f —=~r ex p{~ \ (*- V-) ^ ~H X ~ H)1 cfri ^ 2 A*! 

d (v 27 i) g { j 



where. 



x = 









g 2 = det[2], 



and D = +x£+x£ <, P 2 }, 



. . . . dp 

By the implicit function theorem, > 0 implies that the 



equation P(P;|i,2) = CONSTANT defines a differentiable function 



2 



such that P(R(\i , 2) ; \i , 2) = CONSTANT. If the CONSTANT is 
set to , then i?(jj,,S) is the SEP function. Let X 1 ,X 2 , -,X n be 

a random sample of independent and identically distributed 
random vectors with distribution N(n,2) . Then the maximum 

likelihood estimators of p and 2 are ji = — ^2 X i and 

n i=l 

* « 

2 = — . When ji and 2 replace n and 2 in the 

n i = 1 

equation i^i?;|i,2) = then J?(ji,2) gives SEP, an estimator of 
SEP. 

The second-order approximation for the variance of SEP is 
given by: 

o 2 sBp = (D^SEP) T P(p^SEP) 

+ \ (D 2 y,SEP \ T (F^P) {D^SEP) 

+ ±{ lD^SEP \ T Ej ( 1 ) 

where 2“ = (o 11 , o 12 , o 13 , o 22 , a 23 , o 33 ) r , D^ ^SEP is the derivative of 
SEP with respect to n,2 u viewed as a column vector, D^ U SEP is 



3 



the usual second derivative matrix (or Hessian) of SEP with 



respect to |j.,2 u , 



P = 



n 

P T 
V mS" 




is the variance-covariance 



matrix (1, 2“, ® is the tensor product operation and the 

underline notation is used to represent the column vector 
formed by stringing out the elements of the matrix by rows. 



For example, 



if A = 



1 2 
3 4 



then A = [1 2 3 4] T . 



Note that (1 is 



distributed N( u,— 2) so P„ 

n 



— . It can also be seen thatnS 

n 



has a Wishart distribution with parameter 2 , so that 

a 

= ~ (2<S>2) . The complete derivation of equation (1), based 
on Taylor Series expansion, is given in appendix A. 
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Since (i and 2 are independent [Ref. 4:p. 102] it follows 




Derivation of the first order terms in equation (2) has 
been accomplished by S.D. Hill at JHU-APL [Ref. 5], In order 
to solve for the second order terms it will be necessary to 



find expressions for 



d 2 SEP 

an 2 



a 2 sep , a 2 sep 

and 

as u2 a^as 0 



i.e. it is 



necessary to find D * 2U SEP , the second derivative matrix of SEP 

with respect to its arguments. Chapter II will describe the 
derivation of the second order variance estimate and Chapter 
III will discuss implementation and testing of the algorithm. 
Conclusions will be presented in Chapter IV. 
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II. DERIVATION OF THE SECOND ORDER ESTIMATE 



A. FIRST ORDER TERMS 

We begin by describing the derivation of the first order terms 
given in equation (2) . If we define: 



F[x x ,x 2 ,x 3 ;\i,'Z ir ) = exp J - (x- ji) *2 _1 (x- p) J 

where 



a 11 


o 12 


o 13 


a 12 


= o 21 


a 21 


a 22 


a 23 


, o 13 


tl 

Q 

U) 

H* 


a 31 


a 32 


a 33 


a 23 


= o 32 



(x-h) t 2* 1 (jc-^)= (x 1 -^ 1 ) 2 a 11 +2(x 1 -n 1 )(x 2 -^ 2 )a 21 +2(x 1 -^ 1 )(x 3 -n 3 )o 31 
+ ( x 2 - H 2 ) 2 ° 22 +2 ( x 2 - H 2 )(*3 "^ 3)° 23 + (* 3 - H 3 ) 2 ° 33 



and 

g(p , 0, <J>; p , S u ) = F{psin0cos<{>, psin0sin<|>, pcos0; p, E u ) . 



Then, with a shift to spherical coordinates, we see that: 

P(i?;p,E u ) = — j— f 2 * f*f R p 2 sin6g(p , 0, <t>; p , E u )dpd0d<t> . 

( v /2i) 3 g Jo JoJo ( 3 
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We wish to solve equation ( 3 ) for each of the partials 



dSEP 

dx i 



where X — (Hi/ H2' ^3' °i2' °i3' °22' °23' °33) • From 



1 

previous section, J?(|x,S) = — is the SEP function, so we 



dR 



to solve for the partials . 

0 X ,• 



Taking partials of equation ( 3 ) we see that. 



ffmM. = - — L 



P\R) 



Iff dX, dXldX 2 dX3 



dx 



i (v/27i) 3 gr J 



and 



P'(R) = — El — f 2 *f 1 ’sindg(R, 6, 4>)d0c7<}> . 

(v/ 27 c) 3 g J ° J ° 



We 



i=l, 2,3 



can solve the right hand side of equation 
by using the relationships 

dg __ dg dg_ dg dg __ dg 

dx x 8|i 2 dx 2 8jj. 3 dx 3 



(4) 



the 



need 



(4) 



for 



7 



so we see that 




D 1 



If we apply Green's Theorem to this equation and then make 
a shift to spherical coordinates we get: 

-fffj&dx^dx, = R 2 J' 2 *fJcos4>sin 2 Qg(R,d,<l))ddd<p 

D 1 



-J J J-^-dx 1 dx 2 dx 3 = R 2 f 2n f n sin4>sin 2 dg(R / 0, <J>)d0d<J> 

D 2 0 0 



~III < ^ x i < ^ x 2 < ^ x 3 = R 2 J 2 "J 1 'cosQsinQg(R / Q,<t>)ddd(t > 

D 3 0 0 

= •^-J , o 2 ’ t J , Jsin20g(i?,e / <j))d0d4) 



To solve the right hand side of equation (4) for i=4 , . . . , 9 
we first note that: 



&g_ 

d\i 1 d\i 2 



da ij 

dg_ 

da ij 



i=J 






8 



so we get 



dg 



d °ij 



_l_d__dg_ i= • 
2 dxj 0 ^ 

i + i 

dxj a^ 



which is used to evaluate the integral of 



_dg_ 

ax. 



as follows; 



~fff = ^2 J 2, 'J 1t [ all (P s i n 0 cos 4 ,- K l i) + al2 (Psin0sin4)-|i 2 ) 



+a 13 (pcos0-n 3 )]cos<j)sin 2 0g(i? / 0, 4>)d0d<J) 



-fff •J^dx 1 dx 2 dx 3 = ■R 2 J' 2 ”J , '[o 11 (psin0cos<t)-n 1 )+o 12 (psin0sin4)-|i 2 ) 

D 5 0 0 

+ a 13 (pcos0-p 3 )]sin<t)sin 2 0g(i? / 0, <J>)d0dc{) 



- j J / dx 1 dx 2 dx 3 = j "j' 7 '[o 11 (psin0cos<t)-n 1 )+o 12 (psin0sin4)-p 2 ) 

D 6 0 0 

+a 13 (pcos0-|i 3 )]sin20g(i? / 0, 4>)d0d<t> 

-fff-lfdx^dx, = -^-J 2 ”J' , '[a 21 (psin0cos<{>-n 1 ) + a 22 (psin0sin<j)-|i 2 ) 



+ a 23 (pcos0-n 3 )]sin4>sin 2 0g(i? / 0, 4>)d0d(J) 



"// f c * x id x 2 c t x 3 = J' 2n j'"[o 21 (psin0cos<})-n 1 )+a 22 (psin0sin4)-p 2 ) 

D 8 0 0 

+o 23 (pcos0-|i 3 )]sin20g(i? < 0, 4>)d0d4> 
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^ X 2,t y^ 7l [ a3:L (P s i n ® cos< t )- M' 1 ) + o 32 (p s ±n 0 sin < t > -ti 2 ) 



D 9 

+a 33 (pcos0-jx 3 )]sin209(i? / 0 , <J>)d 0 c?<t> 



If we define the functions: 

cc i.j.k.i( R i t 1 ' S ") = f*cos i (jQ)cos k (l<b)g(R,Q,$)ddd<b 

CS iijikl (R; |i,S“) = j** j* co s i (jB)sin k (l<\>) g(R, 6 ,4>)d&d<\> 

sc i,j,k,i ( R :\ l » s “) = J' 2n J'Jsin 1 (j6)cos k (l<l))g(R, 0 , <j>)d0d<t> 

5S i jk j(i?; p, S") = f 2 ” r”sin i (j0)sin Jc (i<J))sf(i?, 0, <|))d0d<{) 

• J • • Jo Jo 



Extending the work of Blischke and Halpin, with the aid of 
Green's Theorem, we can see that: 



asffifrx, S u ) 



^ 2 , 1 , 1,1 *^ 2 , 1 , 1,1 

•^ < -l ,1,0,0 ^^1 ,1,0,0 



*^1 , 2 , 0,0 
2 ‘^l , 1 , 0,0 



i r 



■*I(/!(m,S u );|i,I u ) 



(5) 
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£££ 2 ,1,1,1 

> 3‘S , s 2 , 1,1,1 



dSEP(\i,l, u ) 

as u 



2 sc. 



1,1, 0,0 



dPi 

d^- 1 , 2 , 0,0 

a^! 

^^ 2 , 1 , 1,1 

dp 2 

d-^-i ,2,0,0 

3P 2 

d-5Ci,2, 0,0 

2^3 



( 6 ) 



l(ii(u,s“);„.r"). 



Applying the partials given in equation (6) to the 



partials - J j J dx 1 dx 2 dx 2 given above we see that: 



dSEP 
ao 11 



1 

2 W 5 



[4o 11 w 9 +2o 12 w 11 +o 13 (w 1 -w 2 )j 

-w 8 (o 11 |x 1 +o 12 |i 2 +o 13 |i 3 )} 



dSEP 



do 



12 



“j-f [ 2 o 11 w 11 + 4 o 12 w 12 +° 13 (w 3 -w 4 )] 
-Wio(o 11 Hi+o 12 ji 2 + o 13 p 3 )} 



dSEP 
ao 13 



^{^[° 11 ( V l' V 2) + ° 12 ( lV, 3-^4) + 0 13 (W 7 + ^5)] 
-w 6 (o 11 iii+o 12 l i 2 -t-o 13 ii 3 )} 



11 



dSEP 



da 



22 



= 2^{f[ 2 ° 12lV, ll +4022,V, 12 + 023 ( V 3 _ ^)] 

-V 1 o(o 12 m + 0 22 ^2 +<j23 » i 3)} 



dSEP 

da 23 



_i_{|[o 12 (w 1 -v 2 )+ o 22 (v 3 -v 4 )+o 23 (v 7 + l / 5 )] 

-v 6 (o 12 m + o 22 Ji 2 + o 23 Ji 3 )} 



dSEP 

ao 33 



_^||[0 1 3( Vl -V 2 ) + 0 23 (V 3 -^) + 0 33 (W 7 +^ 5 )] 

-V 6 ( 0 l 3 J A l + ° 23 J i 2 + 033 l i 3 )} 



where &{i?;p,E u ) is given by: 






W 1 


= 


cc 

^ 1 , 1 , 1,1 


W 2 


= 


^-" 1 , 3 , 1,1 


W 2 


= 


^ 1 , 1 , 1,1 


W 4 


= 


C 5 1 , 3 , 1,1 


W 5 


= 


^^1. 1 , 0,0 


W 6 


= 


SC - i '2' o,o 


w 7 


= 


‘^■' 1 , 3 , 0,0 


W S 


= 


^ 2 . 1 , 1,1 


w 9 


= 


"^ 3 , 1 , 2,1 


W 10 


= 


*^ 2 , 1 , 1,1 


W 11 


= 


^ 3 , 1 , 1, 2 


$ 

to 


= 


■^ 3 , 1 , 2,1 



where each of the functions is evaluated at i?(p. , S u ), p, , S u . 

This completes the derivation of the first order terms 
required in the computation of the variance estimate. 
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B. CONCEPT FOR DERIVATION OF THE SECOND ORDER TERMS 

We first observe that equations (5) and (6) enable one to 

express D * EU SEP as a composition of three functions 

T 1 , T 2 and T z , each of which is a fairly simple function of its 

arguments . 

Explicitly, 

D^^SEF{\ i,E u ) = r 3 oT 2 oT>,S u ) 



where, 



r>,s u ) 



E u 

5EP(n,S u ). 



T 2 (\i,-L u ,R) 



E u 

R 

HR;\i,2 u ). 



13 



2 U 

R 



w o 



w c 



w. 



10 



W, c 



^6 

2w c 



* (4o 1 V 9 +2o 12 w 11 +o 13 (w 1 -w 2 ))-^-(o 1 Vi+o 12 n 2 +o 13 fi 3 ) 



8w c 



2 Wr 



(2o 11 w 11 + 4o 12 w 12 + o 13 (w 3 -w 4 ))--^(o 11 p ;L + 0 12 ii 2 + 0 13 ii 3 ) 



4w c 






W c 



-^-(o 11 (v 1 -w 2 )+a 12 (v 3 -w r 4 )+o 13 (w’ 7 + w 5 ))--^-(a 11 p. 1 +o 12 |i 2 +o 13 n 3 ) 



i? 



w. 



g v '( 2a12 V ll +4a22v i2 + ° 23 ( lvr 3- W 4)) - ‘2^‘( 0l2 ^l +022 | 1 2 +a23 ^ 3 ) 



4w, 



(a 12 (v 1 -v 2 )+a 22 (v 3 -w 4 ) + a 23 (v 7 +w 5 ))--^-^(a 12 ji 1 +a 22 |i 2 +o 23 jA 3 ) 
^-(a 13 (v 1 -w 2 )+o 23 (v 3 -w 4 )+a 33 (w 7 +v 5 ))--^-(a 13 p. 1 +a 23 |i 2 +a 33 n 3 ) 



Since ^ u SEP{\i , 2 U ) = , 2 U ) we may apply the chain 

rule to obtain the second derivative as a matrix product: 



Dl ^SEP{\i , 2 U ) = DT 2 (T 2 or^ll , 2 u )) -DT 2 (T^\l , 2 u )) -DT^il , 2 u ) . ( 7 ) 



The simplicity of the functions T lt T 2 and r 3 means that 
the derivatives DT lt DT 2 and DT 2 are easily obtained and then 
the rather large matrix products involved in DT 2 DT 2 DT^ can be 
left to the computer leading to a rather painless computation 
of Dl yu SEP. 

H/ 2j 
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DT X (the total derivative of the function T x ) is the 10x9 



matrix consisting of the 9x9 identity matrix (from 



c*H,2 u ) ' 



with the tenth row consisting of the partials of SEP with 
respect to mu and sigma. 



DT X 



1 

0 

0 

0 

0 

0 

0 

0 

0 



0 

1 

0 

0 

0 

0 

0 

0 

0 



0 

0 

1 

0 

0 

0 

0 

0 

0 



0 

0 

0 

1 

0 

0 

0 

0 

0 



0 

0 

0 

0 

1 

0 

0 

0 

0 



0 

0 

0 

0 

0 

1 

0 

0 

0 



dSEP dSEP dSEP dSEP dSEP dSEP 



6pi d\i 2 d\i 2 



da 



li 



3a 



12 



3o 



13 



0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
10 0 
0 10 
0 0 1 
dSEP dSEP dSEP 

do 22 do 23 d°33 



DT 2 is the 22x10 matrix consisting of the 10x10 identity 
matrix (from ^ >_Bl \ with rows 11 to 2 2 consisting of the 

12x10 sub-matrix of partials of W with respect to (i, E u and 
R. 
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1 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


i 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1 


dw 1 


dw x 


aw t 


dw ± 


Bw x 


Bw 1 


aw x 


dw 1 


a^ 


dw x 


dpi 


dp 2 


ap 3 


do ll 


a °12 


a° i3 


da 22 


da 23 


da 33 


dR 


9W 12 


SW 12 


dW !2 


dw i2 


dw i2 


dw 12 


dWi2 


dw !2 


dw !2 


dw !2 


api 


dp 2 


ap 3 


da 11 


d °12 


do 13 


d °22 


d °23 


d °33 


dR 



DT 3 is the 9x22 matrix consisting of the partials of 
u i (i=l,... / 9) with respect to p, E u , R and W. 



■ Bu , 


Bu. 


Bu. 


Bu. 


ap (lx3) 


-(1x6) 

as u 




Sw^ xl2) 


Bu 9 


3u 9 


Bug 


BUg 




—(1x6) 

as u 


SR^ X1) 


au< lxl2 \ 



Thus, in order to obtain the second degree variance 
estimate given in equation (2) we must solve for the values of 
all elements of DT lt DT 2 and DT 3 . Once these values are 

computed, the solution is found by the matrix multiplication 
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shown in equation (7) , yielding the 9x9 matrix of second 
partials of SEP with respect to mu and sigma. 

The values of the DT X elements are computed from equations 

(5) and (6) . The determination of DT 2 requires the evaluation 
of and 4^ which is fully described in appendix B. 

Finally, to solve for the elements of DT 3 it is necessary to 

solve for each of the partials given above. The evaluation of 
these partials is given in appendix C. 
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III. IMPLEMENTATION AND TESTING 



A. PROGRAMMING 

The calculation of o^ p the second order approximation for 

variance of SEP (when jj. and 2 are known) has been implemented 

in the FORTRAN program SEPCMP which is provided in appendix D. 

The following steps are followed to obtain the approximation: 

Input the distribution mean (n) and variance- 
covariance matrix elements (2 U ). 

- Compute SEP using p and 2 . 

Compute the values of each element of the matrices 
DTI, DT2 and DT3 and then form the matrix product 
comprising the second derivative of SEP with 
respect to and 2 U . The first derivative of SEP 
with respect to n and 2 U is then given by row 10 
of the matrix DTI. 

- Calculate the variance-covariance matrices of (1 
(P^) and of 2 (P L ) from 2. 

Form the SEP variance approximation using equation 

( 2 ) . 

Computation of SEP given jj. and 2 U is based on the PL-1 
program RAP provided by JHU-APL. RAP is a general program 
which provides the radius of coverage for any probability in 
the interval (0,1) for either 2 or 3 dimensions. Those 
portions of the program applicable to this problem, ie. 
probability = 0.5 and 3 dimensions, have been used in the 
FORTRAN subroutine FINDSEP to obtain SEP. Derivation of the 
equation used to compute SEP is given by L. S. Simpkins [Ref. 
6] and is implemented in the subroutine FINDSEP using Gaussian 
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Quadrature (via the IMSL subroutine QTWODQ) to evaluate the 
required double integral. 

Gaussian Quadrature is also used to evaluate each of the 
double integrals in the W, CC, CS, SC and SS functions, which 
are required in the evaluation of the matrices DTI, DT2 and 
DT3 . 



B. TESTING 

The purpose of deriving the second order approximation of 
the variance of SEP is to yield more accurate confidence 
intervals for SEP. The approximation of these confidence 
intervals is derived from the fact that since p and 2 are 
maximum likelihood estimators of p. and E and thus SEP is 
asymptotically Normal with mean SEP. Thus, for example, an 
approximate 95% confidence interval for SEP is given by 

SEP± 1.96 sJvar(SEP ) . 

In order to test the affect of using the second order 

approximation for the variance of SEP, tests of confidence 

interval coverage were conducted as follows: 

Sample 30 random vectors from 7^n,E) and compute 
(1 and 2 for the sample. 

- Calculate S&P from ft and 2 . 

Calculate d 2 s g p by replacing £ and 2 for p. andE 
in both the first order and second order 
approximations for the variance of SEP. 

Compute the large sample approximation 95% 

confidence intervals for SEP using SEP± 1 . 96<Jo 2 s g p . 
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Repeat 1000 times to calculate the percentage of 
times that confidence intervals based on each 
approximation cover the true SEP. 

For each input variance-covariance matrix, S , 
repeat the test at different magnitudes of the 
mean. 

Trivariate normal random vectors are obtained using the 
IMSL subroutines DCHFAC and DRNMVN. DCHFAC performs Cholesky 
factorization of the distribution SIGMA which is then used as 
an input to DRNMVN which returns the desired random vector. 
The estimates of P and P z are found by replacing p and 

S with (1 and 2 in the formulations of P^ and given in 

Chapter I. The results of this test are summarized in Table 



1 . 





TABLE 1 


. 95 PC'] 


C CONFIDENCE INTERVAL C 


COVERAGE 


1 


SIGMA 




MU 


FIRST ORDER 


SECOND ORDER 










Cl 


Cl 


1.0 


0.8 


0.9 


0 


94 . 6 


95.2 


0.8 


1.0 


0.8 


0 






0.9 


0.8 


1.0 


0 












1.0 


92.8 


93.3 








1.2 












0.8 






100 


-160 


270 


0 


93 • 3 


94.2 


-160 


400 


-480 


0 






270 


-480 


900 


0 












100 


95.6 


95.8 








250 












500 






10,000 


3,000 


2,000 


0 


93.3 


93.8 


3,000 


10,000 


3,000 


0 






2,000 


3 , 000 


10,000 


0 












8,000 


93 . 5 


94.0 








10,000 












12,000 
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These results indicate that there is little accuracy to be 
gained by using the second order variance approximation. 
There is, however, a situation in which the second order 
approximation becomes important. This case has been described 
by K. V. Kitzman [Ref. 3] and is paraphrased here. 

Weapons system targeting can be measured from separate 
subsystems (such as navigation, guidance, postboost, etc.). 
Suppose that there are n independent subsystems and let 

X ilt ..., 2) represent the errors for the subsystem 

i = 1, . . . ,n, then the impact errors for the weapons system are 

n HcL 

Yj = 22 X U j=l, . . . ,m and - N(\i Y =n\i, E r =nE) . 

i= l 

If the Y's are measured directly and used to estimate^ 

and (to get the system SEP) then P u = = — S and 

m m 

p z Y = (S®S). If on the other hand the individual 

n 

X's are available from the subsystems then where 

i= 1 



ana 



Then P u = — S and 
* r m 
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p z Y ~ The significant factor is that the estimate of 

P Sy differs from the first estimate by a factor of n while the 



estimate of P remains unchanged. Thus as the number of 

subsystems becomes larger and the number of observations 
becomes smaller the size of P s can become very small in 

relation to P^, . As the mean approaches 0 it follows that 



dSEP 

dp 



approaches 0 as well and thus the contribution of P^ 



to 



the estimate can become insignificant even though it is large 
in relation to the P 2 term. 

This affect has been examined as follows: 

Divide ^ and S by n, the number of iid subsystems 
being simulated. 

- Draw m random samples from the distribution given 
by this new p. and S for each of the n subsystems. 

- Calculate p and 2 for each of the n subsystems. 

- Sum the n estimates to get a total p and 2 , the 
estimate for mean and variance in Y. 

Compute confidence intervals for SEP using both 
the first and second order approximations with 

P = -2 and P s = — (2®E). 
m ** m v ' 



This test was conducted with n = 4 and m = 5. The results 
of this test are given in Table 2. As can be seen, the effect 
is greatest when the mean is 0, and decreases as the mean 
becomes larger, which is what would be expected. 
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TABLE 2. 95 PCT Cl COVERAGE, MULTIPLE SUBSYSTEMS 



SIGMA 


MU 


FIRST ORDER 
Cl 


SECOND ORDER j 
Cl 


100 -160 270 


0 


92.4 


97 . 6 


1 -160 400 -480 


0 






270 -480 900 


0 








5 


91.4 


94.8 




10 








20 








25 


89.4 


92.2 




50 








80 








50 


92 . 6 


94 . 4 




150 









200 








To test the sensitivity of the second order variance 

approximation procedure to the assumption that the detonations 

are normally distributed, tests were conducted in which random 

vectors were drawn from a mixed normal distribution. The test 

is conducted in the following manner: 

Prior to each sample obtain random variable p from 
U(0,1) distribution. 

- If p is less than .98 then draw a sample from 
and proceed. 

if p is greater than .98 then draw a sample from 
W[|i, 10*11) and proceed. 

Continue until 30 samples have been drawn, then 
compute all values and proceed as above. 

This procedure simulates a distribution with the same mean 
but larger tails than the normal distribution. The confidence 
interval coverages are then compared to those obtained in the 
base case from the same input distribution. The results of 
this test are summarized in Table 3. 
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TABLE 3. 95 PCT Cl COVERAGE, MIXED DISTRIBUTION 



1 - 

1 

i 


SIGMA 




MU 


FIRST ORDER 


SECOND ORDER 


i 








Cl 


Cl 


i 

1.0 


0.8 


0.9 


0 


86.6 


87.6 


0.8 


1.0 


0.8 


0 






0.9 


0.8 


1.0 


0 












1.0 


91.2 


91.8 








1.2 












0.8 






100 


-160 


270 


0 


88.2 


88.4 


-160 


400 


-480 


0 






270 


-480 


900 


0 












100 


91.5 


92.4 








250 












500 






10,000 


3,000 


2,000 


0 


79.3 


80.4 


3,000 


10,000 


3,000 


0 






2 , 000 


3,000 


10,000 


0 








8,000 

10,000 

12,000 


85.6 


86.2 
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IV. CONCLUSIONS 



Examination of the results in Table 1 shows that the 
second order approximation procedure provides a slight 
improvement in the 95% confidence interval coverage. The real 
utility of this method is brought out by the results of the 
multiple subsystem test shown in Table 2, where confidence 
interval coverage improves by as much as 5%. Since this case 
represents the actual situation in missile testing (i.e. very 
small mean and independent system observations) it is 
reasonable to state that use of this procedure will improve 
the estimate of the variance in SEP. Furthermore it seems 
highly unlikely that there would be any advantage to be gained 
by pursuing a higher order estimate of the variance of SEP. 

The mixed distribution test shows that the second order 
approximation is only slightly more robust than the first 
order approximation if the distribution is actually from a 
mixed and not a true normal. It also indicates that if the 
variance is extremely large then the affect of the mixing is 
greater than for the smaller values and that only a few 
extreme values are required to greatly reduce the confidence 
interval for the estimate. 
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APPENDIX A 



In order to derive the equation for the second order 
approximation for variance of SEP (equation (1)) we start by 
noting that the second order Taylor Series Expression for the 
SEP function is given by (ignoring the remainder term): 

SEP(X) = SEP^X)+[D ilJ1 SEP(X)] T AX + ^AX‘ I [D^' L SEP{X)]^X 

where X, D^^SEP, and D^-^SEP are defined in the text, X is the 

vector consisting of the maximum likelihood estimators for 
each element of X and AX = X-X . Our goal then is to find an 

approximation for the variance of SEP(X). Because SEP(X) is 
asymptotically unbiased and the mean square error 

£{(SEF{X)-SEP(X)) 2 ] = Var(SEF(X))+(E{SEF(X)]-SEP(X)) 2 , 
we can get this approximation by looking at E[SEP(X) -SEP(X) ] 2 . 
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Applying the tensor identity ABC = A®C T i2 to the last term 
and then transposing we can write the second order Taylor 



equation as: 



SEP(%)-SEP(\) ~ p (1<E 5SP(>.)] T Al + ^(AX 7 ®AA. :r )[D^_ s SE^ : )] 
= p,, v ,S£flfX)] r A X + -| [Ai 2 . s SB^X)] r (A X® A X) 



Now we note that since £(X) = X it follows that E(AX) = 0 . 
If we square each side and take the expected value we see 
that: 

BpBR^-SBi^X)] 2 = p„ fS »,S£’fl[X)] r £[AA,AA, J ]p |lfl ;u5Efl[X)] 

+ j | p 2 yu 5BP(X)] t £{(A A.® A X)(A X® A X) ^P 2 ^BBPfX) ] 

In order to evaluate £[(AX®AX)(AX®AX) r ) we note that 

var(A X® A X) = 2 (SR) (P®P) ( SR) T [Ref. 7:p. 44] where the matrix 

SR has the following qualities [Ref. 7: pp. 32-38]: 

SR is symmetric 
- (SR) (SR) = SR 

2(5R)(U®^) = U&X + X&U 
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for example, if u and v are of dimension 2 then 



SR 



1 0 

0 — 

2 

0 — 

2 

0 0 



0 0 




0 1 . 



Then since, 

£[(AA®AA)(AA®AA) 1 ] = var(A A® A A) +£(A X® A A)[£(A A.® A A .)] T 



and 



AAA A T = A A.® A A. 



if we define P to be £[AA,AA. r ) = var(AA) we see that 



var( A A.® A X) +E( A A® A A)[£i;A A® A A )] T 



2 SR(PS>P)(SR) t +ee 1 
2 (SR) T (Pg>P)(SR) +E El . 
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Since D^^SEP is a symmetric matrix it follows that 
(SR) Dy v,SEP = D 2 ?SEP and we see that: 

* p^ E Stfi^)] r ip KfS SOT(A.)] 

+ [d 2 , s S£HX)] T P2)J ^D U 2 . sSgfffl) ] 

+ 1 T EE ^Du-nSEP(X) ^ , 

which in turn is approximately equal to the variance of 
SEP ( X ) . 

Finally, since S is a 3x3 symmetric matrix, we need only 
describe the unique terms (which we have denoted by E u ) to 
arrive at the second order equation for SEP given by equation 
( 1 ) • 
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APPENDIX B 



A. APPROACH 

It is seen that each of the 12 W i functions consists of 

a double integral of cos and sin functions multiplied by the 
function g(.Rsin0cos4>, Rsin0sin4>, Rcos0) . Each of the required 
derivatives carries through the integrals, so that we need to 
evaluate each in respect to the function g. 

We will first evaluate each of the partials of g, and then 
carry these through to the functions W i . 

Recall , 

= exp|--i[(i?sin0cos4>-n 1 ) 2 o 11 

+ 2(i?sin0cos<|>-|i 1 )(i?sin0sin<|>-ji 2 )a 12 
+ 2(f?sin0cos<t>-|i 1 )(i?cos0-p 3 )o 13 
+ (i?sin0sin4>-p, 2 ) 2 o 22 
+ 2(Rsin0sin<t>-n 2 )(i?cos0-ii 3 )o 23 
+ (i?cos0-^ 3 ) 2 o 33 ]} . 
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If we define the vectors A lf A 2 , A 3 and A 4 by 



Ax 



sin 2 0cos 2 <j) 

2sin 2 0cos<J>sin<t> 

2sin0cos0cos<}> 

sin 2 0sin 2 <}> 

2sin0cos0sin<J> 

cos 2 0 



^2 



sinOcos^m 

sin0sin<t>n 1 +sin0cos4)n 2 
cos0|i 1 +sin0cos4>n 3 
sin0sin<J>n 2 
cos 2 +sin0sin<}>|j.3 
cos0n 3 



= [m>i 2 Hi^2 2 ^i^3 ^2 2h 2 h 3 H3) t , 



^4 



•RsinOcos^-i^ 
i?sin0sin<t>~n 2 . 
i?COS0-(i 3 
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Then treating R as an independent parameter and taking the 
partial with respect to R we have: 

= -gr-{a 11 (i?sin 2 6cos 2 <|)-p. 1 sin0cos<|)) 

+cr 12 (2.Rsin 2 0cos<}>sin4>-n 1 sin0sin<t>-p 2 sin0cos4>) 

+a 13 (2i?sin0cos0cos<t>-n 1 cos0-n 3 sin0cos(j)) 

+a 22 (Rsin 2 0sin 2 <t>-n 2 sin0sin<f>) 

+a 23 (2Rsin0sin4>cos0-n 2 cos0-n 3 sin0sin<t>) 

+a 33 (i?cos 2 0-p. 3 cos0)} 



Similarly, we derive: 



dg 

d\i ± 



= g r *{a 11 (i?sin0cos<t>-n 1 )+a 12 (i?sin0sin4)-n 2 )+a 13 (i?cos0-n 3 )} 
= g-jfa 11 a 12 a 13 ]^ 4 } 



= sr-{CT 12 (Rsin0cos<t>-n 1 )+a 22 (Rsin0sin<t)-|i 2 )+a 23 (Rcos0-|i. 3 )} 
= g-{[a 12 a 22 a 23 ]^ 4 } 



= s r -{a 13 (Rsin0cos4>-ji 1 )+o 23 (Rsin0sin<t>-p 2 )+a 33 (Rcos0-n 3 )} 
= sr*{[a 13 o 23 a 33 ]A 4 } 
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In order to derive the partials with respect to sigma we 

is ij 

will first qenerate the partials and the partials of 

°°kl 

with respect to a kl . These values are given below. 





nil 


ell 


nil 


a 22 


a 

0011 


-o 11 © 11 


- 0 12 a ll 


-a 13 o 11 


° 33 -a 22 © 11 

g 2 


a 

0o i2 


-2a 11 a 12 


”° 33 -2a 12 o 12 

Q 2 


° 23 -2o 13 o 12 

g 2 


-2o 22 a 12 


a 

0Ol3 


-2a 11 a 13 


-^-2a 12 a 13 

q 2 


■° 22 -2o 13 o 13 

g 2 


_20 13_2 0 22 0 

g 2 


a 

d°22 


°33 _ a ll a 22 

Q 2 


-a 12 a 22 


~° 13 -o 13 a 22 

g 2 


-o 22 o 22 


a 

0o 23 


~ 2 °23 - 2 o u 0 23 


-^i 3 -20 12 0 23 
g 2 


° 12 -2o 13 o 23 

g 2 


-2a 22 o 23 


a 

do 33 


° 22 -o 11 a 33 

g 2 


"° 12 -o 12 o 33 

g 2 


-o 13 o 33 


0ll -a 22 o 33 

g 2 





nil 


nil 


XJCi 


a 


_<J 23 -23-11 


H22 -a 33 o 11 


-o 11 


aaii 


g 2 


g 2 


2g 


a 


° 13 -?n 23 o 12 


_2a i2 _ 0a 33 a 12 


-o 12 


^°12 


g 2 


g 2 


g 


a 


°12 -2q23„13 


-2o 33 o 13 


-o 13 


0013 


— £a\J U 

g 2 




g 


a 


_ a 23 a 22 


Cl1 -o 33 o 22 


-a 22 


d°22 


O U 


g 2 


2g 


a 


"°1 1 -->-23-23 


-2o 33 o 23 


-o 23 


5 °23 


g 2 




g 


a 


— O 23 O 33 


-a 33 o 33 


-o 33 


a°33 






2g 
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We can now use these values to find the desired partials 



dg - 



da 



11 



= ^^[(iesinecos^-jx^o 11 ) 2 

+2(i?sin0cos<t)-|a 1 )(i?sin0sin<t>-n 2 )a 12 o 11 
+2(i?sin0cos4>-(i 1 )(i?cos0-|i3)o 13 o :l1 

+(i?sin 0 sin<})-n 2 ) 2 ^a 22 o 11 --^^ j 
+2(J?sin0sin4>-n 2 )(i?cos0-ia. 3 )^o 



23 0 ll + ^23 





) 


L33 0 H_ <T 22 > | 


O 11 


[ Q 2 l 


2 1 



gr-|-| [a- 1 l R- 2 (S -i) % - 2 Ra 1 1 (S - 1 ) U A 2 + a 11 (S - 1 ) U A 2 ] 

[R^- 2RA 2 +A 3 ]\ 



0 0 0 



33 23 



g 2 g 2 g 2 



= g* 



0 H( S - 1 ) U+ o o 0 



~°33 °23 "°22 



,11 



[R 2 A r 2 RA^A 2 ]-^- { 



Following this same methodology it is easy to see that: 




34 



3a. 



|o 23 (E _1 ) u +( ° 23 ° 13 


1 

Q 

> K 

o 


°n 


o) 


fi? 2 A 1 -2i?A 2 +A 3 l-a 23 i 


N 

b 1 

oq 

N 

fcF 


2g 2 


2g 2 


1 





Ja- = g-\i 

aa 33 y [2 



0 3 3 ( S -i)U + | _^22 _£l2 0 °11 o o 

g 2 g 2 g 2 



[R 2 A 1 -2RA 2 +A 3 \— 
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We are now able to calculate each of the partials of thefts 

by simply computing the effects of the four A vectors on the 
original W functions and applying the formulas derived on the 
previous pages. We will carry through the complete derivation 
for W 1 and then show the solutions for the other cases. Since 

the vector A 3 does not involve trigonometric functions, it 
will not be necessary to compute any partials for it. 
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